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I. 


INTRODUCTION 


A  problem  of  considerable  current  interest  is  the  modeling  of  laser 
propagation  through  inhomogeneous  media.  Of  particular  interest  is  the 
case  where  the  index-of-refraction  inhomogeneities  are  induced  by  high- 
energy  beams  intez acting  with  the  medium.  The  self-consistent  modeling 
of  the  problem  requires  the  modeling  of  the  hydrodynamics  of  this  inter¬ 
action,  the  effect  of  the  hydrodynamics  on  the  index  field,  and  the 
influence  of  this  induced  field  on  the  propagation.  Such  models  have 
been  derived  from  first  principles  by  a  number  of  authors.1"6*  The  effect 
of  the  phenomenon  has  been  examined  by  use  of  analytical  approximations7*8 
and  by  numerical  solution  employing  a  number  of  algorithms.1"^1** ***  The 
most  successful  algorithm  now  used  by  several  of  the  large-scale  computer 
codes  is  the  so-called  Implicit-Galerkin-Spline  (IGS)  technique  devised 
by  Herrman  and  Bradley.9  Even  more  important  than  the  numerical  algo¬ 
rithm,  however,  is  the  form  in  which  the  model  is  preprocessed  for  compu¬ 
tation.  Experience  indicates  that  focused  beams  of  high  Fresnel  number 
require  the  model's  being  specified  in  an  adaptive  coordinate  system; 
otherwise,  as  the  marching  computations  approach  the  focus,  spatial  reso¬ 
lution  is  insrfficient  for  any  numerical  scheme  utilizing  a  reasonable 
number  of  grid  points.  Furthermore,  the  focusing  (defocusing)  strategy 
must  allow  for  the  fact  that  the  beam  will  be  distorted  and  spread  to  an 
extent  usually  lying  somewhere  in  expanse  between  the  exit-plane  spot 
size  and  the  so-called  diffraction-limited  spot  size  at  the  focus.  For 
high  Fresnel  numbers  the  transformations  to  adaptive  coordinates  must  be 
accompanied  by  an  appropriate  phase  transformation;  otherwise,  phase 
oscillations  cannot  be  sampled  adequately.  Transformations  previously 
used  to  overcome  this  difficulty  either  partially  removed  the  initial 
phas*'  oscillation  and  allowed  arbitrary  coordinate  focus  (defocus) 3. 5 
or  removed  all  the  initial  phase  but  were  limited  to  j^focusing  strategy 
following  the  behavior  of  a  beam  focused  in  vacuum. 

The  generalized  transformations  derived  herein  contain  the  best 
features  of  both  these  cited  methods  and  are  shown  to  be  a  convenient 
point  of  departure  for  solution  by  a  class  of  numerical  methods.  Two 
computer  codes  utilizing  these  transformations  have  been  developed  and 
compared,  one  using  a  modified  form  of  the  Herrmann- Brad ley  algorithm 
and  the  other  using  a  Fast  Fourier  Transform  (FFT)  approach  described 
herein,  /he  formulation  cited  above  is  shown  to  lead  to  an  FFT  solution 
that  does  not  require  a  Nyquist  accuracy  criterion,  allowing  the  numeri¬ 
cal  procedure  to  march  the  solution  forward  in  a  more  economical  fashion. 
Numerical  experiments  indicate  that  the  IGS  method  is  slightly  faster  on 
a  per  step  basis  when  identical  grid  spacing  is  used.  The  FFT  approach, 
however,  yields  greater  accuracy  for  identical  spacing  of  mesh  in  the 


* 

References  'ire  listed  on  page  41. 

**Thie  list  is  representative  but  by  no  means  complete. 

***See  vote  after  Reference  9. 


7 


Preceding  page  blank 


) 


transverse  plane  and,  because  of  larger  allowable  propagation  steps,  it 
is  overall  faster.  This  factor,  in  conjunction  with  the  greater  ease 
of  implementing  an  FFT  algorithm,  arj  the  requirement  for  less  core 
storage  (when  symmetry  is  not  assumed),  makes  the  method  superior  to 
algorithms  currently  used. 

The  algorithms  described  herein  concern  themselves  primarily  with 
the  solution  of  the  propagation  equation.  The  class  of  problems 
addressed  is  represented  by  hydrodynamic  equations  which  can  be  solved 
essentially  in  closed  form  within  each  computing  step.  This  is  the 
type  of  problem  to  which  most  attention  has  been  directed  in  the 
literature. 

Stimulus  for  the  FFT  portion  of  this  report  arose  from  a  three-way 
consultation  between  the  author,  Dr.  J.  Wallace  of  Far  Field,  Inc., 
Sudbury,  Massachusetts,  and  Dr.  J.  Lilly  of  the  US  Army  Missile  Command. 
A  related  FFT  technique  for  the  wave  equation  and  hydrodynamics  for 
repetitively  pulsed  lasers  is  described  in  an  Army  Missile  Command 
Report 22  which  is  in  preparation.  The  nature  of  this  collaboration  is 
described  in  Section  V  and  the  Acknowledgment. 


II.  MATHEMATICAL  MODEL 

We  consider  the  model  applicable  to  continuous-wave  operation  of  a 
laser  in  steady-state  conditions.  We  assume  that  the  laser  is  propagating 
through  an  atmosphere  with  a  time -independent  cross  wind  V(£)  and  is 
being  rotated  (slewed)  at  an  angular  rate  Q  through  an  axis  passing 
through  the  exit  aperture.  We  assume  that  the  hydrodynamics  can  be 
represented  by  the  steady-state  heat  equation  with  the  laser  acting  as  a 
source  term,  that  the  wind  and  slewing  produce  forced  convection,  and 
that  conductivity  is  negligible.  Under  these  assumptions  the  governing 
model  is  given  by* 

2ik  3W/3C  ♦  (32W/3C2  +  32W/3n2)  +  k2(r.2  -  nQ2)W  =  0,  (2-1) 

n2  -  nQ2  =  (nQ2  -  l)Ap/p  =  -(n02  -  1)AT/TQ,  (2-2) 

Vp(y  +  ai  e  1  lwlZ  *  (2-3) 

fc'U.n.O)  --  (P/7ir2)FU,n)exp[i$  -  i(k/R)(£2  +  n2)/2].  (2-4) 


Equation  (2-1)  is  the  paraxial  approximation  to  the  wave  equation. 
Equation  (2-2)  is  obtained  from  the  lorentz-Lorenz  law  expressing  the 
change  in  the  index  of  refraction  as  a  function  of  the  change  in  tempera¬ 
ture  or  air  density  where  small  deviations  from  ambient  are  assumed. 
Equation  (2-3)  is  the  expression  of  the  heat  balance  for  the  conditions 
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stated  above.  Equation  (2-4)  expresses  the  amplitude  and  phase  of  the 
launched  beam.  In  (2-4)  the  two  quantities  in  the  complex  exponential 
represent  the  phase  arising  from  the  laser  cavity  and  the  focusing 
optics,  respectively.  It  should  be  noted  that  in  Equation  (2-3)  we  have 
neglected  the  so-called  kinetic  cooling  phenomenon,^3  since  this  effect 
is  insignificant  for  ground-based  lasers.  This  effect  can  nevertheless 
be  brought  into  this  model  without  change  in  the  basic  nature  of  the 
algorithm,  since  it  merely  affects  a  quadrature  formula. *>3>3 

We  convert  the  model  to  a  form  containing  all  variables  in  dimen¬ 
sionless  form  by  letting 


u  = 

W/CP/irr2)55, 

(2-5) 

X  = 

C/r,  Y  =  n/r. 

(2-6) 

z  = 

Ul, 

(2-7) 

a  = 

a.L, 

(2-8) 

6  = 

kr2/L, 

(2-9) 

Y  = 

(nQ2  -1)  P/(T0P0cpVXr), 

(2-10) 

P  = 

I  kL(n2  -  no2), 

(2-11) 

O)  = 

AL/V, 

(2-12) 

f(X,Y)  = 

F(Xr,Yr), 

(2-13) 

*(X,Y)  = 

$(Xr,Yr). 

(2-14) 

Here  we  have  assumed  that  r  is  a  convenient  scale  length  in  the 
transverse  (£,n)  plane.  For  a  „aussian  beam  a  convenient  choice  for  r 
is  the  e”1  folding  radius.  For  focused  beams  the  characteristic  length, 
L,  in  the  ;  direction,  is  conveniently  chosen  to  be  the  focal  cistance  R. 
In  order  to  make  one  set  of  equations  valid  for  focused  beams  and  colli¬ 
mated  beams  (R=«°) ,  L  is  defined  by 


L  = 


R,  focused  beams 
kr2,  collimated  beams. 


(.2-15) 


For  this  choice  of  dimensionless  variables  the  model  takes  the  form 


2iB3U/3Z  +  (32U/3X2  +  32U/3Y2)  +  26pU  =  0, 

f* 

H(X,Y,Z)  =  -aye'aZ(l  +  <oZ )  l  J  |U(x‘  ,Y,Z)  1 2dX*  , 


(2-16) 

(2-17) 
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U(X,Y,0)  =  f(X,Y)  exp(i*  -  ig^X2  +  Y2)/2],  (2-18) 


where 

g  = 

(kr2/R,  focused  beams 

(2-19) 

(  1  ,  collimated  beams, 

and 

gl  = 

|  g ,  focused  beams 

(2-20) 

{  0,  collimated  beams. 

The  solution  to  the  model  expresses  the  dimensionless  irradiance 
lur  as  a  function  of  the  four  dimensionless  parameters  a,  g,y,u)  and 
X,Y,Z.  If  we  seek  the  maximum  irradiance  or  a  characteristic  average 
irradiance  in  the  focal  plane,  then  the  solution  is  seen  to  be  dependent 
only  on  a,g,y  and  <*>. 

Equation  (2-16)  is  of  parabolic  type,  and  the  algorithms  success¬ 
fully  utilized  for  solution  have  similarities  to  those  employed  for 
numerical  solution  of  the  time-dependent  heat  equation.  The  Z  variable, 
representing  the  propagation  direction,  is  analogous  to  the  time  variable 
in  heat  transfer.  The  various  methods  used  thus  far  "march"  the  compu¬ 
tations  forward  in  Z  as  is  done  in  the  classical  explicit  and  implicit 
difference  schemes.  Equation  (2-18)  is  used  for  specifying  the  initial 
plane  of  data.  At  each  step  the  index  distribution  is  recalculated, 
utilizing  the  most  current  values  for  U  where  AT,  as  will  be  shown,  is 
evaluated  by  numerical  quadrature.  Despite  these  similarities  to  the 
heat  transfer  problem,  the  algorithms  referenced  above  have,  of  neces¬ 
sity,  been  uniquely  tailored  to  the  propagation  problem  and  are  somewhat 
more  complicated  than  the  techniques  used  for  solving  the  classical  heat 
equation.  Factors  of  computing  efficiency  of  10-100  can  be  gained  by 
transforming  the  above  model  into  a  form  more  suitable  for  computation. 
These  transformations  are  essential  regardless  of  the  numerical  algorithm 
ultimately  used. 


III.  TRANSFORMATIONS 


Adaptive  coordinates  are  generally  required  for  focused  beams  of 
high  Fresnel  number  but  have  also  been  found  necessary  for  collimated 
beams  when  significant  beam  spreading  occurs.  The  transformations  uti¬ 
lized  for  adaptive  coordinates  and  phase  removal  are  as  follows: 


x  =  X/a(Z) ,  y  =  Y/a(Z), 


(3-1) 


U  =  [B/a(Z) ] 


exp[y  ig  (X2  +  Y2) (da/dZ) /a (Z) ] , 


z 


(3a2)_1dZ  . 


(3-2) 

(3-3) 
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The  application  of  these  transformations  to  (2-16)  results  in 


2i  3B/3z  +  (32B/3x2  +  32B/3y2)  +  g(x,y,Z)B  =  0  (3-4) 

where 

g(x,y,Z)  =  2Ba2y  -  B2(x2  y2)a332a/3Z2.  (3-5) 

The  choice  of  the  term  (3a/3Z)/a  in  the  phase  transformation  (3-2)  makes 
possible  the  maintenance  of  a  simple  structure  for  the  transformed  propa¬ 
gation  equation  for  any  subsequent  choice  of  scale  length  a(Z).  All 
first-order  derivatives  B  ,  B  that  would  otherwise  appear  in  the  trans¬ 
formed  equation  have  beenxel¥minated  by  this  particular  choice  of  vari¬ 
ables.  The  transformations  (3-1),  (3-2),  and  (3-3)  are  a  generalization 
of  transformations  used  by  Wallace, 2  Aicken  et  al.,3  and  Herrmann  and 
Bradley. ^  By  appropriate  choice  of  a,  these  cited  transformations  result 
from  the  above  generalization.  In  a  later  section  we  will  show  how  cer¬ 
tain  advantages  peculiar  to  each  of  the  cited  system  of  transformations 
can  be  obtained  by  appropriate  choice  of  a(Z). 

The  solution  to  (3-4)  is  highly  oscillatory  in  z,  portending  diffi¬ 
culty  in  numerical  computation  due  to  the  Nyquist  criterion.  The  oscil¬ 
lations  arising  from  the  term  containing  g  are  removed  by  use  of  the 
Herrmann- Brad  ley9  transformation 

D  =  B  exp(-ir/2)  (3-6) 

where  z 

r  =  f  g(x,y,Z)dz.  (3-7) 

zo 

This  leads  to  the  equation 

2i  3D/3z  =  -exp(-  y  ir)(32/3x2  +  32/3y2)[D  exp(-  y  ir)].  (3-8) 

If  we  envision  the  approximation  of  (3-8)  by  finite  differences,  where 
all  difference  operators  are  centered  at  z  =  z0  =  zn  +  Azn/2,  then  all 

terms  containing  T  become  unity.  The  resulting  difference  equation  is 
then  identical  to  the  difference  equation  arising  for  a  collimated  beam 
in  vacuum.  The  difference  approximation  to  Equation  (3-8)  is  then  indis¬ 
tinguishable  from  the  approximation  to 


2i  30/ 3 z  =  -(32/3x2  +  32/3y2)D.  (3-9) 


In  sequence  with  the  transformations  leading  to  it,  Equation  (3-9)  is 
the  point  of  departure  for  numerical  solution  using  both  the  IGS  and  FFT 
mevhods. 
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IV.  IMPLICIT-GALERKIN-SPLINE  METHOD 


The  IGS  algorithm  is  designed  to  advance  the  solution  of  Equation 
(3-9)  from  the  plane  zn  to  the  plane  zn+1  =  zn  +  Az.  The  algorithm,  in 
conjunction  with  the  initial  wave  function  specification,  is  applied 
repetitively,  along  with  phase  initializations  at  each  step,  to  march 
the  solution  to  the  point  in  z  of  interest.  Since  we  have  assumed  the 
wind  to  be  along  the  x-axis,  we  have  a  line  of  symmetry  along  y=0.  We 
use  this  symmetry  to  reduce  the  calculations  in  half.  In  the  rectangu¬ 
lar  propagating  tube 


o  <.y<\. 

0  <  J  <  :M, 


(4-1) 


we  superimpose  the  three-dimensional  grid 


Xj  —  “X|^  (j-1)  Ax,  j  —  1,2, ..  .J 


Xk  =  (k-l) Ay, 


z„^,  =  z„  +  Az. 

n+l  n  n’ 


k  =  1,2, ...K 
n  =  1,2,. ..M. 


(4-2) 


On  this  grid  we  approximate  Equation  (3-9)  by  a  Crank-Nicholson  type 
stencil  in  z  and  the  resulting  difference  equation  by  a  fractional 
step-alternating  direction  approximation.  This  leads  to  the  two 
equations : 


(1  -  j  iAz32/3x2)Dn+Jj  =  (1  +  ~  iAz32/3x2)Dn, 


(4-3) 


(1  -  ~  iAz32/9y2)Dn+1  =  (1  +  j  iAz32/3y2)Dn+?5, 


(4-4) 


where  the  superscripts  denote  the  wave  function  at  the  n,  n+k,  and  n+l 
z  levels,  respectively. 
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The  linear  spline  approximation  is  given  by 


D  ~  D  = 


where 


1= 1  m=l 


^(x)wm(y) 


(4-5) 
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0, 

X  f  X^J 

or  x  >  xuv 

m  ,j 

(x-x^1)(x^-xM) 

*l-i  -  x  5  *v 

4*1,  J 

(4-6) 

<x*i 

-  x)  (xitl  ■  xt) 

xzix<  xt^. 

4*i,J 

»!(*) 

1  ,  2 
*  exp  2  <X1 

-  x2),  -«  <  X  <  X 

(4-7) 

*j(x) 

1  /  2 
*  exp  j  ( Xj 

-  X2),  Xj  <  X  <  ». 

(4-8) 

A  similar  definition  applies  to  “n(y)* 

24 

Applying  the  Galerkin  method  to  (4-3),  one  obtains 

00 

/ w.  (x)  ((1  -  j  iAz32/3x2)ffn+,S  -  (1  ♦  ~  iAz32/3x2)Fn]dx  «  0,  (4-9) 

•  00 

with  a  similar  expression  resulting  from  (4-4).  After  performing  the 
integration,  one  obtains  two  sets  of  coupled  linear  equations  given  by 


(l*iai)CD;.iik.D;.1>k).(4-2i.l)D"k 


(4-10) 


where 


.  D^)  .  (4*2i.,)Dj;‘ 


aj  =  3/2  Az/h2, 
h  =  Ax  «  Ay. 


(4-11) 


(4-12) 

(4-13) 


It  should  be  noted  that  the  effect  of  using  the  Galerkin-Spline  procedure 
as  opposed  to  the  more  conventional  central  difference  method  has  been  to 
yield  a  finite  difference  stencil  of  similar  basic  structure  but  differing 
coefficients.  The  central  difference  approach  would  yield  (4-10)  and 
(4-11)  with  the  coefficients  (1  ±  iaj)  being  changed  to  (iiaj)  and 
(4  t  2iaj)  being  changed  to  (6  ±  2iaj).  The  use  of  the  Galerkin-Spline 
is  reported ly^  more  accurate. 
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Equations  (4-10)  and  (4-11)  are  of  tri -diagonal  fora  and  are  solved 
by  the  recursive  relations  obtained  by  an  adaptation  of  the  Gauss  elimi¬ 
nation  aethod.  The  complete  set  of  equations  for  both  the  x  and  y  direc¬ 
tions  are: 


x  Direction 


A  *  C  *  1  -  iaj 

(4-14) 

Bj  *  4  ♦  2ia1#  j  ■  2,3,...J-1 

(4-15) 

Bj  «  Bj  «  2  ♦  3(1  -  o.sx^2)/xjq 
♦  a j (1  *  O.SXj^h  ♦  0.25h(l  -  0. SXj^"2)^”1) i , 

(4-16) 

V  “  *  i.,)(D;.1>k  ♦  Dj,lk)  *  t4  -  2i»,)Dj>k, 

(4-17) 

E,  •  -C/Bj , 

(4-18) 

P,  ■  G,/*,. 

(4-19) 

Ej  ■  -w  Ej-.  *  v“- 

(4-20) 

Fj  *  <G}  •  A  V.)(A  Ej.,  *  »,)'*.  J-2.3.-J 

(4-21) 

DJ*1  k  ’  »• 

(4-22) 

“"Ik  *  EjD$,k  *  Fj-  j  ■  . 

(4-23) 

The  sequence  of  relations.  Equations  (4-14)  to  (4-23),  when  performed 
for  k  =  1,2,...K,  comprises  the  first  half  of  the  fractional  step  method. 
It  should  be  noted  that  all  quantities  are  complex  and  that  complex 
arithmetic  is  implied  throughout. 

y  Direction 

The  algorithm  for  the  y  direction  is  similar  to  that  for 
tion  with  changes  arising  due  to  the  assumed  line  of  symmetry 

the  x  direc 

» 

A  *  C  *  1  -  ialf 

(4-24) 

=  4  +  2ialf  k  =  1,2, ...K-l 

(4-25) 

Bk  *  Bj  Eq.  (4-16) 

(4-26) 
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I  A  *  rZ7**7*  *gg 


1 


\ 

l 

’ 

s 


} 

t 


i 

t 

i 

[; 


*  -  *•  x*'.***-,  * *•  .v.-^sh—  -  »  >. «■■•#-  vv*-  ,  •»*. 

> 

I 

} 

*> 


* 

•)l 


(4-27) 

i 

i 

t 

El  -  -2C/BJ  , 

(4-28) 

X 

i 

5 

Fi  «  Gi/Bi, 

(4-29) 

* 

Ek  «  -C(AE|c_1  ♦  Bk)'\  k  *  2,3, ..  .K 

(4-30) 

Fk  -  (Gk  -  A  FkM)(A  Ek_j  ♦  8k) ” 1 ,  k  -  2,3,... K 

(4-31) 

dM.  i  *  "• 

(4-32) 

■  EkC.i  *  Fk-  k  * 

(4-33) 

Upon  completion  of  this  second  half  of  the  alternating  direction  method, 
the  array  d"  has  been  replaced  with  D1?*1.  When  combined  with  the  phase 

}lV  J  A- 

initialization.  Section  VI,  this  operation  advances  the  solution  one 
step.  A  FORTRAN  IV  subroutine  for  this  algorithm  is  listed  in  Appendix  B. 


V.  FAST  FOURIER  TRANSFORM  METHOD 

Fast  Fourier  Transform  algorithms  have  been  employed  for  solving 
boundary  value  problems  by  a  number  of  investigators. 25-27  These  efforts 
involved  the  use  of  the  FFT  to  reduce  the  dimensionality  of  certain  dif¬ 
ference  equations,  the  resulting  numerical  methods  requiring  a  mixed 
procedure  entailing  aspects  of  the  more  traditional  implicit  difference 
schemes  as  in  Section  IV.  An  FFT  algorithm  for  solving  the  atmospheric 
propagation  problem  for  collimated  beams  was  proposed  by  Wallace,19  and 
led  to  the  effort  being  reported  here  and  in  the  MICOM  report22  pre¬ 
viously  mentioned.  An  FFT  algorithm  for  propagation  within  a  laser 
cavity  is  described  by  Canavan  and  Phelps.2*1*21  The  authors  report  on 
one  method  applicable  to  collimated  beams  which  utilize  the  Helmholtz 
equation  as  a  point  of  departure.  A  second  technique  is  reported  which 
is  applicable  for  diverging  beams  and  which  makes  use  of  the  paraxial 
approximation.  This  latter  form  contains  a  singularity  at  z=0  and  does 
not  seem  to  be  generally  applicable.  The  algorithm  proposed  by  Wallace19 
is  designed  for  the  solution  of  the  paraxial  equation  for  collimated 
beams  and  is  based  on  a  forward  difference  approximation  in  the  Z  direc¬ 
tion.  This  algorithm  would  march  one  step  of  the  solution  to  Equation 
(2-16)  by  the  three-stage  operation 
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U(p,q,Z)  «  FFT[U(X.Y,Z)eiewAZ]  , 


(5-1) 


U(p,q,Z*AZ)  *  U(p,q,Z) (1  ♦  iAZ(p2  ♦  q2)^]'1,  (5-2) 

U(X,Y,Z  ♦  AZ)  »  FFTj[U(p,q,Z*  AZ)],  (5-3) 

where  the  symbolism  FFT  and  FFT.  indicate  operation  on  a  two-dimensional 
array  of  complex  data  by  the  forward  and  inverse  FFT  algorithm.  This 
algorithm,  while  being  competitive  with  existing  methods,  suffers  the 
disadvantage  of  requiring  a  Nyquist-type  accuracy  criterion.  For 
example,  when  u*0  (vacuum  propagation),  the  exact  FFT  transform  solution 
to  (2-16)  requires  that  in  (5-2)  IT  be  given  by 

U(p,q,Z  ♦  AZ)  =  U(p,q,Z)e‘lAZ(p2  +  q2)/2.  (5-4) 

Since  this  solution  is  cyclic,  the  Nyquist  criterion  indicates  that  a 
minimum  of  six  samples  per  cycle  is  needed  to  integrate  a  harmonic  of 
a  particular  frequency.  When  Equation  (5-2)  is  used,  the  Nyquist  cri¬ 
terion,  when  applied  to  the  highest  frequency  (it/ Ax  *  n/Ay),  leads  to 

AZ  <  2Ax2/3*.  (5-5) 

For  the  IGS  method  Herrmann  and  Bradley  suggest  the  forward  marching  be 
governed  by  the  step-size  criterion,  AZ  «  Ax2,  which  is  approximately 
one  sample  per  cycle  for  the  Nyquist  frequency.  As  they  show,  however, 
this  criterion  does  lead  to  sizable  but  tolerable  phase  errors  at  the 
high  frequencies.  Hence,  in  (5-5)  and  (5-2)  one  could  similarly  sacri¬ 
fice  the  accuracy  of  the  high-frequency  components  by  relaxing  (5-5). 
From  this  the  conclusion  can  be  drawn  that  the  algorithm  given  by  (5-1) 
to  (5-3)  requires  a  step-size  criterion  that  is  similar  to  the  IGS 
method  but  differs  by  a  factor  near  unity  to  be  determined  by  numerical 
experiment.  The  two  algorithms  thus  differ  in  their  cost  according  to 
their  relative  cost  in  marching  the  solution  one  step  forward,  other 
factors  being  assumed  equal.  As  will  be  detailed  in  a  later  section, 
the  single-step  IGS  algorithm  is  about  15  percent  faster  than  the  FFT 
as  executed  on  BRL's  BRLESC  II  computer.  In  contrast,  it  will  be  shown 
by  example  that  the  FFT  algorithm  is  more  accurate  for  a  specified  grid 
size  than  is  the  IGS  method. 

By  analogy  with  the  algorithm  to  be  described  in  a  later  section, 
the  above  algorithm  can  be  made  more  attractive  by  replacing  (5-2)  with 
(5-4).  This  can  be  justified  by  virtue  of  tb*  approximation 


e-iAZ(p2+q2)/2 


[1  +  iAZ(p2  +  q2)/2] 


(5-6) 
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or  by  assuming  constancy  of  u  within  a  step  Az.  By  this  artifice  the 
Nyquist  requirement  is  eliminated;  however,  the  algorithm  can  be  shown 
to  propagate  th*  non-linear  effect  to  an  accuracy  that  is  first-order 
accurate  in  AZ. 


In  order  to  obtain  a  aore  accurate  algorithm,  to  eliminate  the  need 
for  the  Nyquist  criterion,  and,  aore  iaportantly,  to  handle  focused  as 
well  as  colliaated  beaas,  we  derive  an  FFT  algorithm  utilizing  the 
sequence  of  transformations  in  Section  III. 

The  FFT  algoritha  is  also  designed  to  advance  the  solution  to 
Equation  (3-9)  in  stages  froa  plane  zn  to  zn<M,  each  stage  preceded  by 
a  phase  initialization.  For  reasons  peculiar  to  FFT  computer  subroutines 
we  do  not  use  the  line  of  syaaetry  but  work  on  the  rectangle 


-^  <  x  <  ^  , 

-V  *  *  v 


(5-7) 


with  :he  grid  points  as  defined  by  (4-2).  On  this  region,  at  the  arbi¬ 
trary  plane  z,  the  conplex  wave  function  D  is  assumed  to  have  a  Fourier 
series  approxiaation 


where 


D(x,y,z) 


V'T  djk(x,y,z), 


(5-8) 


djk(x,y,z)  «  Cjk(z)  exp(ipjX  ♦  iqky).  (5-9) 


By  assuming  that  each  component  of  the  series  satisfies  (3-9)  in  the 
interval  (z_,zn+1),  we  guarantee  that  (5-8)  also  is  a  solution.  Hence, 
the  amplitude  of  each  harmonic  for  advancing  from  zn  to  zn+1  is  giver,  by 

Cjk(W  *  Cjk(2n}  exP[’  7  iA*(Pj2Mk2)].  (5-10) 


where,  because  of  the  transformations  leading  to  (3-9),  we  solve  the  one- 
step  propagator  "exactly. "  From  the  theory  of  finite  Fourier  series  for 
discrete  data  it  can  be  shown  that 


where 


(JO 


■±t 

1*1  m«l 


o(x£,yn,zn)  exp(-ip^Xj  -  iqmyk),  (5-11) 


!2n(£  -  l)/(JAx) ,  l  *  1,2,..  .J/2 
2n(£  -  1  -  J)/(JAx) ,  l  =  J/2+1, . . .J 


(5-12) 
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(5-13) 


(  2*(m-l)/(KAy),  m  ■  1,2, . .  .K/2 
(  2w(m-  1  -  K)/(KAy),  ■  «  K/2+1, . . . K. 


The  inverse  finite  Fourier  transform  at  level  (n*l)  is  given  by 


D(x»y*Vi) 


J  K 

xEZVVi}  cxp(ipjx  +  iqKy) 

7*1  k-1 


(5-14) 


The  operations  defined  by  the  suns  in(S-ll)  and  (5-14)  can  be 
related  to  the  operations  perforned  by  the  so-called  Fast  Fourier  Trans¬ 
forms.  See  Appendix  A.  In  (5-11)  and  (5-14)  we  have  adopted  a  conven¬ 
tion  which  "folds"  the  negative  frequencies  to  the  second  half  of  the 
series.  This  form  is  compatible  with  the  convention  used  in  most  com¬ 
puter  subroutines.  See,  for  example,  Brenner. 28 


VI.  PHASE  INITIALIZATION 


The  employment  of  the  phase  transformation  (3-6),  (3-7)  is  equiva¬ 
lent  to  propagating  one  z  increment  holding  the  phase  distribution  fixed 
at  what  is  anticipated  will  be  the  distribution  at  the  middle  of  the 
step.  This  "anticipation"  is  implemented  by  assuming  that  the  function  g 
can  be  separated  into  similarity  type  products  of  functions  depending  on 
z  and  on  x,y.  In  this  manner  the  z  dependence  can  be  accounted  for  by 
quadrature,  and  within  the  half  step  the  x,y  dependence  is  assumed  to  be 
similar  to  the  behavior  at  the  beginning  of  the  step.  Once  the  step  has 
been  taken,  by  application  of  the  basic  algorithm,  the  phase  of  the  wave 
function  is  now  "off"  by  one-half  step.  The  phase  is  updated  by  again 
assuming  separability  of  g  and  similarity,  this  time,  however,  using  the 
x,y  dependence  of  g  calculated  at  the  end  of  the  step.  This  phase  up¬ 
dating  can  be  combined  with  the  phase  initialization  for  the  next  half 
step  and  done  in  one  operation.  If  the  scaling  parameter  a  is  completely 
arbitrary--that  is,  its  slope  is  discontinuous  at  each  step--then  an 
additional  phase  increment  must  be  included.  These  operations  are 
defined  in  the  equations  that  follow.  We  let 


where 


and 


Z)  = 

f  [fjCZMx.y^)  ♦  f2 

(Z)(x2  ♦  y2)]dz. 

(6-1) 

zn+*s 

*l(Z) 

=  2a6Ya[Z(z)]e’aZ  [l  + 

uZ(z)]"1, 

(6-2) 

=  -BV^a/'OZ2 , 

(6-3) 

u(x,y,z)  =  -J  |D(x,y,z)| 

2dx. 

(6-4) 

—  00 
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1  ?SS}^ .. ..  ^  v.'  «*>»>« 


The  phase  updating  is  accomplished  by  multiplying  the  "old"  wave  func¬ 
tion  by  a  complex  phase  increment  and  storing  the  "new"  wave  function 
over  the  old;  that  is 


>  *  Dold(x’y’2n>e’U** 

(6-5) 

where 

^-Ty(x,y,zn)(En*Fn) 

+  £  c^hvvv-v* 

fz  n 

(6-6) 

E"’ j 

i 

f  fl(Z)dz, 

2n-JsA2n-i 
(Z  n 

(6-7) 

rn’J 

J 

'  f 1 (Z)dz , 

V^n 
f  2n 

(6-8) 

C„-J 

f  f2(Z)dz, 

zn‘^zn-i 

fz  n 

(6-9) 

"n'J 

f  f2(Z)dz, 

ln+Jjain 

(6-10) 

I  «  - 
n 

Ba(2n)3a/3Z  .  . 

n 

(6-11) 

J  * 
n 

Ba(Zn)3a/3Z|ZsZn+  . 

(6-12) 

In  the  computer  code  the  four  integrals  (6-7)  to  (6-10)  are  done  by 
trapezoidal  quadrature.  When  the  scaling  length  a  has  a  continuous 
slope  in  going  from  step  to  step,  the  sura  (In  +  Jn)  vanishes. 


VII.  CHOICE  OF  SCALING  LENGTH  FOR  ADAPTIVE  COORDINATES 

The  scale  length  a(Z)  used  in  defining  the  adaptive  coordinates  is 
present  in  the  phase  transformation  (3-2).  The  initial  wave  function 
after  transformation  is  given  by 

i(4>+d>1) 

B(x,y,0)  =  af (ax,ay)e  (7-1) 

where 

♦  jfx.y.O)— i  (x2  ♦  y2)(63a/3Z|z=0  +  fjj).  (7-2) 
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Unless  a(Z)  is  properly  chosen,  +1  can  become  quite  large  away  from  the 
center  of  the  beam,  leading  to  rapid  oscillations  in  the  initial  data. 
Such  oscillations,  to  be  sampled  adequately,  would  require  too  fine  a 
mesh  size  for  practical  computation.  The  manner  in  which  this  diffi¬ 
culty  is  avoided  is  described  below  for  focused  beams  and  collimated 
beams. 

Focused  Beams 


For  focused  beams  Bj  -  6  in  Equation  (7-2).  Hence,  if  3a/3Z  =  -1 
at  Z  =  0,  will  be  identically  0  as  desired.  The  choice  of 

a(Z)  *  [1-Z)2  ♦  (Z/B)2]*5  (7-3) 

by  Herrmann  and  Bradley®  and  by  Ulrich*-*  has  this  advantage  but  suffers 
the  disadvantage  of  focusing  as  does  a  beam  in  vacuum,  thus  not  allowing 
enough  grid  space  for  severely  bloomed  beams.  The  choice  of 

a  *  1  -  (1  -  N/B)Z  (7-4) 

is  similar  to  Aitken  et  al.  and  has  the  advantage  of  being  N  times  dif¬ 
fraction  limited  at  the  focus.  Here  N  can  be  chosen  arbitrarily, 
depending  on  the  expected  severity  of  thermal  blooming.  This  choice  has 
the  disadvantage  that  4^  does  not  vanish  completely;  that  is 

^(x,y)  *  N(x2  ♦  y2)/2.  (7-5) 

Hence,  if  the  N  required  to  contain  the  bloomed  beam  is  too  large,  this 
formulation  may  lead  to  a  poor  numerical  solution.  This  is  apparent  in 
Problems  1  and  2  for  an  N  of  5.  The  above  two  candidates  assume  that 
one  formula  governs  the  z  variation  of  a  from  aperture  to  focus.  This 
requires  that  z0  in  Equation  (3-3)  be  zero.  A  generalization  of  this 
choice  is  to  advance  to  the  focus  in  a  manner  where  a(z)  is  piecewise 
linear.  In  this  way,  a  can  have  an  initial  slope  of  -1,  thus  nulling 
(fij,  and  can  be  gradually  varied  at  each  step  so  as  to  have  the  proper 
scaling  according  to  the  demands  of  the  blooming.  Here  we  have 

a(Z)  -  a„  *  bn(Z  -  Zn),  Z„  <  Z  <  Z„M. 

and  2 

2  =  e'7  [an  .  bncz  -  Zn)l'2dZ 

^n 

«V'V‘  -  Ian*  VZ-Vr‘l-  V0 

(Bany‘(Z-Zn).  b„=0. 

An  additional  attractive  candidate  for  a  is 


(7-6) 


(7-7) 
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a  «  I  -  Z  ♦  NZ2/8. 


(7-8) 


It  has  the  desired  properties  of  -1  initial  slope  and  arbitrary  size  at 
the  focus.  Its  aejor  disadvantage  is  that  it  does  not  permit  the  inte¬ 
gral  in  (23)  to  be  done  in  simple  closed  form. 


The  choice 


a  *  ((1  -  Z)2  ♦  (NZ/Bn 


(7-9) 


has  all  the  requisite  characteristics.  This  form  is  seen  to  be  similar 
to  the  Z -dependent  beam  waist  for  a  focused  Gaussian  beam  with  the  scale 
gradually  evolving  to  N  times  larger  than  that  expected  in  vacuum.  While 
this  form  seems  to  be  an  obvious  generalization  of  the  previously  used9 
(7-3),  the  freedom  to  use  the  generalization  is  not  readily  apparent 
until  one  proceeds  to  perform  the  sequence  of  transformations  as  done 
in  Sections  III  and  VII.  For  a  in  (7-9),  Equation  (3-3)  relating  z  and 
Z  becomes 

z  =  N‘1tan"1[NZ/B(l  -  Z)]  ,  (7-10) 


Z  =  8  tan  z/(N  ♦  8  tan  z). 


(7-11) 


The  U.o  equations  (7-10)  and  (7-11)  thus  permit  the  integration  or 
marching  to  proceed  in  z  space  while  simultaneously  computing  such  quan¬ 
tities  as  f'i(Z)  and  f2(Z)  which  are  more  conveniently  expressed  in  Z. 

It  should  be  noted  that  there  are  "tradeoffs"  among  these  various 
choices.  For  example,  when  a(Z)  is  piecewise  linear,  its  second  deriva¬ 
tive  vanishes,  eliminating  one  term  in  the  phase  function  g(z),  but 
introduces  an  additional  term  because  of  the  discontinuity.  In  Equation 
(7-9)  the  second  derivative  of  a,  as  required  in  (3-5),  takes  the  form 


32a/9z2 


(N2/62)/a3. 


(7-12) 


Choice  of  Az 


Marching  or  integrating  forward  in  steps  of  fixed  Az  is  good 
strategy  since,  through  the  transformations,  the  effect  is  to  use  more 
z  planes  in  the  region  where  the  beam  is  tightly  focused  and  fewer  planes 
elsewhere.  It  can  be  shown  that  for  a(Z)  chosen  by  (7-4)  the  focus  is 
reached  when 


z  =  1/N. 


(7-13) 


Hence,  if  we  elect  to  march  to  the  focus  in  M  equal  steps,  Az  is  given 
by 

Az  =  1/MN.  (7-14) 


For  a(Z)  in  (7-9)  the  focal  z  is  equal  to  u/2N  and 

Az  =  tt/2MN. 


(7-15) 
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In  order  to  have  a  flexible  code  the  author's  codes  have  been  designed 
to  allow  optionally  for  fixed  Az  or  fixed  AZ  with  provisions  for  a 
maximum  AZ  when  Az  is  fixed. 

Collimated  Beams 

For  collimated  beams  01  in  Equation  (7-2)  is  zero.  If  the  compu¬ 
tations  are  done  without  an  adaptive  coordinate  system,  3a/ 3Z  *  0  and 
<t>1  is  identically  zero.  If  the  blooming  is  severe  enough  to  spread  the 
beam  beyond  the  original  coordinate  system,  the  linear  form  of  a. 

Equation  (7-4),  is  generally  adequate.  If  the  required  N  is  not  too 
large  (2  or  3),  then  the  oscillations  can  be  adequately  sampled.  A 
seemingly  better  choice  (though  not  tried)  is 

a(Z)  =  1  ♦  dZ2.  (7-16) 

The  term  <f>j  vanishes  for  this  choice  and 

z  =  2  (2d2  +  2dZ2)"1  ♦  (ZdrHan'^z/dV  (7-17) 

This  form  requires  marching  in  steps  of  fixed  Z,  since  inverting  of  (7-14) 
to  reference  z  dependent  quantities  requires  a  numerical  procedure.  The 
parameter  d  is  free  to  be  chosen  according  to  the  thermal  blooming  that 
is  expected. 

VIII.  t  COMPARISON  OF  THE  ACCURACY  AND  EFFICIENCY 
OF  THE  TWO  METHODS 

The  two  algorithms  were  compared  employing  two  different  methods 
for  choosing  the  adaptive  coordinates  and  utilizing  different  step  size 
Az.  Two  sample  problems  were  done,  one  with  a  "small"  Fresnel  number  of 
2.96  and  the  second  with  a  "moderate"  Fresnel  number  of  9.3.  Both 
cases  experience  severe  thermal  blooming  and  tax  the  state  of  the  art 
for  such  computations. 

Problem  1 

The  first  problem  is  the  so-called  NRL  sample  problem.  We  assume 
an  ideal  Gaussian  beam  focused  at  a  range  of  2  km.  The  parameters  are 
given  by 

f(£,n)  =  exp[-  j  (£2  +  n2)/r2] 
ctj  =  0.064/m 
P/M  =  44.843  kw-sec/m 
r  =  10  cm 
R  =  2  km 
Q  =  0.0/sec 
A  =  10.6  x  10*6m 
♦  =  .0 
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The  parameters  associated  with  the  computations  are  as  follows.  For 
both  methods  the  extent  of  the  grid  at  the  exit  plane  was  ±4  beam  radii 
with  symmetry  utilized  for  the  IGS  method.  An  N  of  5  was  used  for  the 
coordinate  system.  For  the  IGS  method  the  number  of  grid  points  was 
65  *  33  and  64  *  64  for  FFT.  The  actual  size  of  the  individual  cells 
was  essentially  identical  but  twice  as  many  cells  were  used  for  the  FFT 
because  of  the  requirements  of  the  subroutine.  With  this  gridding,  one 
step  of  IGS  was  about  15  percent  faster  than  one  step  of  FFT.  By  impli¬ 
cation,  in  problems  not  permitting  the  symmetry  the  FFT  would  be  over 
1.5  times  faster  than  the  IGS.  Alternatively,  special  modification  of 
the  FFT  subroutines  would  gain  the  same  factor. 

For  this  problem  the  number  of  z  steps  to  the  focus  was  varied  to 
see  if  either  method  permitted  more  rapid  advance.  In  Figures  1  through 
4  the  notations  land  II  are  used  to  denote  the  strategy  for  choosing  the 
adaptive  coordinate  system,  the  notation  implying 

I,  a (Z)  determined  by  Equation  (7-4), 

II,  a(Z).  determined  by  Equation  (7-9). 

For  this  problem  we  have  compared  only  the  peak  intensity  at  the  various 
ranges.  The  solid  line  in  Figures  1-4  is  the  result  obtained  by  Hogge5 
for  the  identical  problem  using  an  explicit  difference  scheme  with  a  grid 
of  121  x  121  in  x,y  and  247  z  planes.  In  order  to  compare  with  Hogge's 
result,  we  have  multiplied  the  intensity  by  eaZ  to  conform  with  his  con¬ 
vention  (private  communication). 

In  Figures  1  and  2,  the  IGS  method  is  seen  to  agree  closely  with 
Hogge  for  10  z  steps,  but  for  5  z  steps  and  a(Z)  chosen  by  I  the  IGS 
method  is  off  considerably.  For  the  FFT  method.  Figures  3  and  4,  the 
results  again  agree  very  well;  but,  most  importantly,  the  solution  shows 
almost  no  sensitivity  to  relaxing  Az  to  5  steps  when  choosing  a(Z)  by 
either  I  or  II.  No  attempt  was  made  to  solve  this  problem  with  either 
method  for  less  than  5  steps.  Neglecting  the  solutions  using  I,  the 
results  for  this  problem  are  inconclusive.  However,  we  note  that  we 
neglected  the  usual  step-size  criterion  in  the  IGS  method.  We  would 
therefore  conclude  that  for  this  problem  the  two  methods  are  comparable 
only  if  we  know  a  priori  that  the  step-size  criterion  can  be  ignored. 

In  Problem  2  we  will  show  that  it  cannot  in  general  be  ignored  for  IGS 
but  can  be  for  the  FFT. 

Problem  2 


The  second  problem  has  more  severe  thermal  blooming  than  Problem  1. 
Again  we  consider  an  ideal  Gaussian  beam  focused  at  2  km.  The  parameters 
for  this  problem  are 

dj  =  0.224/m 
P/V  =  50  kw-sec/m 
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a(Z)  -  II 
•  AZ  =  0.1  n/2N 


D.8  1.2  1.6  2  0 

RANGE  (km) 


em  1  -  Effect  of  Changing  az 
GS  Method 


r  *  25//2  cm 
R  «  2  km 
a  •  0.0/sec 
A  »  10.6  x  1C"6W 
N  *  5 

The  problem  was  solved  utilizing  both  methods,  two  strategies  for  the 
adaptive  grid,  and  various  choices  of  Az.  The  summary  of  the  focal  plane 
data  and  the  key  parameters  used  for  each  of  the  nine  cases  are  listed 
in  Table  I.  Figure  5  denotes  the  isoirradiance  contours  at  the  aperture 
plane  for  each  of  the  nine  cases.  The  focal  plane  isoirradiance  contours 
obtained  by  using  a  subroutine  developed  by  Hartwig^S  are  contained  in 
Figures  6-14.  In  these  figures  the  point  of  peak  intensity  is  denoted 
with  an  x,  and  the  ratio  of  focal  plane  to  aperture  plane  peak  intensity 
is  denoted  as  MAX  FLUX.  The  intensity  for  each  contour  line  is  in 
decreasing  increments  of  0.1  times  the  peak  intensity  as  one  moves  away 
from  the  x.  The  dashed  circle  is  the  e_1  contour  for  vacuum  propagation 
and  for  no  thermal  blooming  encompasses  the  first  six  contours.  In  the 
focal  plane  the  extent  of  the  computing  grid  is  the  original  factor  of  4 
beam  radii  (4  times  the  small  circle)  times  the  factor  of  N*5  to  allow 
for  blooming.  The  magnitude  of  the  focal  plane  blooming  is  evident  by 
comparing  the  area  of  the  inner  six  contours  with  the  small  dashed  circle 
in  Figures  6-14. 

Since  the  problem  has  no  independent  solution  by  which  we  can  gauge 
the  accuracy  of  the  numerical  solution,  we  have  taken  the  usual  approach 
of  trying  to  find  a  region  in  mesh  spacing  where  the  solution  is  not 
sensitive.  From  experience,  accurate  solutions  can  be  delineated  by 
examination  of  the  degree  of  smoothness  in  the  isoirradiance  contours. 
Employing  these  two  factors,  we  have  made  the  judgment  that  the  cases 
denoted  with  the  symbol  ++  are  acceptable  and  accurate.  It  should  be 
noted  that,  since  the  peak  intensity  is  a  quantity  determined  at  a  dis¬ 
crete  point,  some  discrepancy  between  methods  is  to  be  expected,  depending 
on  the  coarseness  of  the  grid. 

For  this  problem  it  was  found  that  a  grid  of  65  *  33  was  not  suffi¬ 
cient  to  obtain  an  accurate  solution  by  the  IGS  method.  Figures  6  and  8. 
For  a  grid  of  101  *  51  with  the  IGS  method  and  a(Z)  chosen  by  II 
(Figure  9),  the  solution  at  the  focal  plane  agreed  with  the  FFT  cases 
using  a  grid  of  64  *  64.  For  the  IGS  method  it  was  found  that  using  the 
step-size  criterion  Az  =  Ax2  yielded  an  accurate  solution  (Figure  9), 
but  when  this  restriction  was  ignored  and  10  z  steps  were  used  (Figure  11) 
an  inaccurate  solution  resulted.  For  the  FFT  method  no  appreciable  dif¬ 
ference  resulted  between  10  z  steps  and  5  z  steps.  Figures  13  and  14. 

The  contrasting  behavior  of  the  two  methods  when  Az  is  increased  is 
attributed  to  the  Nyquist-type  accuracy  criterion  as  discussed  in 
Section  V.  Even  though  implicit  schemes  are  generally  considered  to  be 
unconditionally  stable,  because  of  the  cyclic  nature,  the  Nyquist 
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721  ALPHA  «  0.4480 

.00  BETA  -  8.2618 

1.00  GAW1A  -  44.0375 


150- IRRADI ANCE  CONTOURS 


RUN  NO.  719  ALPHA  >  0. 4480 

Z/L2  •  1.00  SETA  •  8-2616 

MAX. FLUX*  1.60  CANNA  -  44. 0375 
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I  SO- I RRADI ANCE  CONTOURS 


DIANCE  CONTOURS 


RUN  NO.  72 1  ALPHA  »  0.4480 

Z/L2  •  1.00  BETA  *•  9.2618 

MAX. FLUX*  0.77  GAMMA  »  44.0375 


I5Q-IRRADIANCE  CONTOURS 


ISO- IRRADI ANCE  CONTOURS 


IRRADIANCE  CONTOURS 


ALPHA  »  0. 4480 

BETA  •  9.26(8 

GAflHA  •  44.0375 


Method  -  Phase  Oscillations  Not  Removed 


RUN  NO.  723  ALPHA  -  0.4460 

Z/L2  •  I . 00  BETA  >  8.2616 

NAX.FLUX-  0.60  GANNA  >  44.0375 


IRRADI ANCE  CONTOURS 


criterion  is  needed  to  sample  the  oscillations  of  the  high  frequencies. 

A  similar  behavior  would  be  expected  for  the  FFT  method  using  (5-2)  when 
Az  is  relaxed.  The  difficulty  does  not  occur  using  (5-10),  since  the 
oscillations  are  followed  "exactly"  for  each  frequency  within  the  step. 
This  does  not  imply  that  the  solution  itself  is  exact,  since  the  phase 
transformations  done  at  each  step  are  obviously  approximate. 


IX.  CONCLUSIONS 

The  numerical  solution  of  the  laser  propagation  problem  can  be 
greatly  enhanced  by  mathematical  transformations  that  are  applied  con¬ 
currently.  Transformations  previously  used  can  be  generalized  to  remove 
phase  oscillations  arising  from  the  focusing  optics  and  still  permit 
beams  blooming  well  beyond  a  one  times  diffraction  limited  coordinate 
system.  The  FFT  and  IGS  approaches  to  the  solution  of  the  propagation 
problem  are  found  to  differ  most  importantly  in  the  fact  that  a  step- 
size  criterion  similar  to  the  Nyquist  criterion  is  needed  for  IGS  and 
can  be  relaxed  for  certain  formulations  of  the  FFT.  Additionally,  the 
FFT  method  yields  greater  accuracy  than  IGS  for  the  same  choice  of 
gridding.  While  single  steps  of  IGS  are  about  15  percent  more  efficient 
than  FFT,  the  above  factors  led  to  a  factor  of  10  in  computing  efficiency 
for  FFT  for  one  sample  problem,  while  for  another  problem  the  computing 
costs  were  about  equal.  Overall,  the  FFT  method  seems  to  be  several 
times  more  efficient  than  the  IGS  method.  Propagation  programs  based  on 
the  FFT  method  are  easier  to  implement  because  of  the  general  availability 
of  FFT  subroutines.  This  latter  advantage  is  eliminated,  however,  when 
an  IGS  one-step  subroutine  is  available.  See  Appendix  B.  The  IGS  method 
requires  the  storage  of  two  complex  arrays  while  FFT  requires  one.  This 
is  a  significant  saving  of  core  storage  for  those  problems  that  do  not 
permit  symmetry. 
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APPENDIX  A 


USE  OF  FAST  FOURIER  TRANSFORM  SUBROUTINES 


The  operation  performed  by  a  so-called  FFT  subroutine  is  the  very 
rapid  computation  of  the  complex  array 

J 


C.  =  J'1^^  exp[vi2ir(£-  l)(j  -  1)/J] 


(A-l) 


where  the  C-,  j  *  1,2,...J,  are  usually  stored  over  the  complex  array  D^. 
For  v  =  -lathis  corresponds  to  the  forward  finite  Fourier  transform 
over  discrete  data.  This  is  equivalent  to  obtaining  the  amplitudes  of 
the  J  Fourier  harmonics  for  approximating  D.  The  approximation  for  D 
is  then  given  by 


DM  *  ^  cj  exp[vipjX] 


(A- 2) 


where  Pj  is  given  by  Equation  (5-12)  of  the  report  and  v  =  1.  If  x  is 
replaced  by  the  discrete  values  over  which  D»  was  defined,  it  car.  be 
shown  that 


=  =  Cj  exp[viPj(£-  1) Ax] 

J  j  =  1 

=  exPlvi27,(j  -1)(^-1)/J]. 


(A-3) 


The  implication  of  this  identity  is  that  the  sums  (A-l)  and  A-3)  can  be 
performed  by  one  algorithm  with  the  appropriate  change  of  v,  as  is  done 
by  the  various  available  FFT  routines.  The  possible  pitfall  in  the 
propagation  problem  is  that,  in  performing  the  integration  of  the  wave 
function  as  approximated  by  the  Fourier  series,  we  must  use  the  first 
sum  in  (A-3)  rather  than  the  second.  Use  of  the  convention  standardized 
by  the  computer  routines  can  obscure  the  fact  that  the  second  sum  of 
(A-3)  is  a  completely  erroneous  representation  of  the  Fourier  series  at 
all  points  in  x  except  the  discrete  points.  The  difficulty  is  avoided 
by  observing  that  all  operations  involving  x  at  other  than  the  discrete 
points  must  treat  the  series  as  having  the  negative  frequencies  "folded" 
as  given  by  Equation  (5-12).  These  considerations  affect  only  Equation 
(5-10)  of  the  algorithm.  Equations  (5-11)  and  (5-14)  involve  situations 
where  the  sums  of  (A-3)  are  interchangeable. 

The  nature  of  the  original  Cooley-Tukey  FFT  algorithm  and  its  many 
variants  is  such  that  the  operation  is  far  more  efficient  for  arrays  of 
length  The  computer  code  was  designed  to  work  with  arrays  of  size 
64  *  64  to  take  advantage  of  this  economy. 
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APPENDIX  B 


I.  LISTING  OF  FORTRAN  IV  SUBROUTINE  FOR  ONE  STEP  IGS  AOVANCE 


8 

c 

c 

c 

c 

8 

c 

c 


1 .00 


.-Go 


>5u 


IANCI IXM«JYM.DR.OIVOHR(OHIVDZOH2,B1RB.B1IB) 
LOl.5U.DK 101 til ). DHR f 101 . 51). OHI (101.51) 

IT )•  JYM-MIN  REPORT) 


OF  REPORT 


ANO  14-511 
OF  DENOMINATOR 
RECOMPUTATION. 


IMPLIES  X  DIRECTION.  1D*2  IMPLIES  Y  DIRECTION 


SUBROUTINE  AOVANCUXMj 
DIMENSION  DRIlOl,  ~ 

DIMENSION  EKRC 

DIMENSION  PSIRl _ _ , _ 

IXM*J(  IN  REPORT).  JVM-M  IN  REPORT) 

DR.DI  TW j  2  DIMENSIONAL  ARRAYS  FOR  STORING  REAL  AND  IMAGINARY 
PARTS  OF  0  RESPECTIVELY.  STORAGE  APPLIES  AT  LEVELS  N  AND  N+l. 
OHR.DHI  SANE  STORAGE  AS  DR  AND  01  APPLYING  TO  »„5Vel  N+l/2 
0Z0H2»0I/H**2 

B1RB.B1I8  REAL  ANO  IMG  PART  OF  Bi.  EOS. 14-16).  14-26) 

EKR.  EK!  REAL  AND  IMG  PART  OF  E  EQS.(4-18)  ANO  (4-30) 

FKR.  FP.I  REAL  ANO  IMG  PART  OF  F  EOS.  (4-21)  AND  (4-3)) 

PS  II.  PS IR  REAL  AND  IMG  PART  OF  RECIPROCAL 
CONTAINED  IN  BOTH  E  AND  F.  STORED  TO  AVOID 
AIR-1. 

AlI«-i.5*0Z0H2 
81R*4. 

B1!*3.*DZ0H2 
ClR*i. 

C^I*-).5*OZOH2 

10*0 
ID* 10+1 

I F  < I0«c0.l)KlMAX*JYM 
IFUD.E0.2)KiMAX*IXM 
IF! I0.£Q.1)KMAX*IXM 
I F ( ID.tW «2)KMAX*JVM 
DO  100,  K1*1,K1MAX 
IF! lO.EQ  *1)J*K1 
I F ( I0.EQ.2)  I*K1 

INITIALIZE  EK  AND  FK  FOR  RECURSION 

GOT  (200,3001,  ID 

TQR--CIR 

TQ  I*-C1I 

TQDR-B.RH 

TQCI*B1IB 

SUBROUTINE  GK  EVALUATES  GJ  ANO  GK,  EOS.  14-17)  ANDI4-27)  OF  REPORT 
CALL  GKIGKR.GKI, l, 1, J,DZ0H2,B1RB,B1IB) 

GOTO  330 
TQR— 2.*C1R 
TQ  I=-2.*C1 I 
TQDR*B1R 
TQCI=BiI 

CALL  GMuKR,  GKI.L,  I, 1,DZ wH2,BlRB,HllB) 

CONTINUE 

DEN=TOOR*TQOR+TOOI*TgDI 
TR*(Tg*ATQDa*Tgi*Tgnii/PEN 
TI*ITQI*TQDR-TOR*TgDI i/DEN 
EKRI1I-TR 

ek  k  i  )>r  i 

FKRI i )*( GKR*TODR+GK I*TQDI )/UEN 
FK  II 1  )*( GK I*TQDR-GKR*TgDI l/UEN 
00  500  K=2, KMAX 
COMPUTE  EK  ANO  FK 
IF  I 1D.EU • 1) I*K 
IF(  ID.LU.2U-K 
IFIKi.NE.UGOTO  375 
BR*P1R 

IF  I K. Eg. K MAX  JBR=H1RB 

IF (K.  EC. KM AX )BI=BiIB 

CENR*AlR*EKRIK-l)-AlI*EKI(K-lMtJK 

DEN  I  * A iR*EK I ( K- 1 ) +A1IPEKR IK-1 ) ♦ B I 

OEN-=DENR*DENR*DEN  I*OEN  I 

PS  IRIK  J-CENR/OEN 

PS  I II K )*-OEN I/DCN 

EK  R(  K  ?=-ClR*PS IRIK  )+ClI*PSI K  K) 

tK  II K )=-C  1R*PS  IKK  )— Cl Z*PSIR<  K ) 
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,75 


">wU 


'J  W 


I  \>\J 
t'OU 


LL  GK( GKR, GK I* ID. I , J,DZ0H2t 81RB 

<=gkk-air*fkr(k-IuaIi*fkiJk-ii 

I  =  GK  I-A1R*FK I( K-i)-AlI*FKRtK“il 


CALL 
TOR 

TQI  . .  . . . 

FKR(M*TUR*PSIR(K)-rUl«PSIIIK) 
FKHK)  =  TUR*PSII(K  UTU^PSIRIK) 
EKR( KMAX )».o 
cK  II  KMAX  )*«<J 


B1RB.B1IB) 


UKR  ». 

UKI  >  =  .  .‘ 

co  :ooj  kk=i,kmax 

K=KMAX*  .-KK 

UKR=  EKRI K )*UKRO-EX l ( K )*UK I0*FKR (K ) 

UK  I  =  tKR(K  )*UK IO*EK I(K )*UKR0+FKI  IK) 

GOT -•(  6  j  !•*  70J ) »  in 

CHR(K,J)=UKR 

CH I ( K » J ) =UK  I 

GOTO  800 

CR(  I.KMUKR 

Cl (  f tK)=UKI 

CONTINUc 

UKR  =UKR 

UK  I  =UK  l 

CONTINUc 

IF(  IO.cO.DGJTu  iOo 

RETURN 

fcNC 


L 

C 


Jw 


0  j 


SUBRjUTI  it  GK(GKk,GKI#IO, If J,0Z:H2,H1RB,B1IB) 

COMM^NIUSE  MAIN) 

SUBROUTINE  EVALUATES  GJ  ANO  GK,  EQS.  (4-17)  AND  (4-27)  uF 

REP  RT  AS  NEEDEC  ev  ADVANC  FOR  IGS 

AlR-i. 

A1I=  ; . ^  *DZ  JH2 
131R  =  4. 

B1  I*-i,*0Z0H2 
Cl R-  i  • 

Cl  1=  1 • 3  *DZ  jH2 
GO  T  j  (  *0  •  .  20.'),  10 
IF (  I.LG.IlULR*. 

IF ( I.Nt.llULR*  OR! I-1,J) 

IF (  I .  EG.  i  )UL  1  =  ., 

IF  (  I.NE.  .  )UL  (-  (HI 
UCR-  UR<  I,J) 

UCI  =  i; I (  ( •  J  > 

1 F  (  I.I.G.  IXM  )URR=.  . 

IF  (  I.N-.  IXM)URR= 

I F  (  I .  FO .  I XM )  UR  I  =  .  .. 

I F  C  I  .Nr .  1XM)URI  =  OH  1  +  1,  J) 
eR*c:.« 

Bh"il 

IF  I  I  .EG.  i.OK.I.,0.  IXM  )BR*B1RB 
IF( I.LJ. l.nu.l.CG. IXM)bI=-BiIB 


GOT.'  i'.  1 

IFU.Ew.  i  )UL  R*DHR  (  IVJ«1I 
IF  I J  .  EG.  .  )UL  I  =0H  1 1  I,J+U 
IFIJ.Nw.  1  )ULR=DHK ( I , J- l ) 
IF(J.N£.i»UL I  =  DH 1 1 I.J-ll 
UCR=OHR( l,J) 

UC I  =  CH I (  I, J) 

IF ( J . bO. JYM IURR- . ' 

IF!  J.Nr:.JYM)URR=f/HR| 
IFIJ.FG.  IYM)URI  =  .C 
IFIJ.NF. JVMIURIsDHIl I*J«ll 


br*m;r 

B I *w i  I 

I F I J.EQ.JYM>BR=«1RH 
IF ( J • EG. J YM ) B I *-H  i  IB 
CONTINUE 

GK R  =  A  1R*ULR- A 1  I*iJL  I*BR*UCR-BI  *UC I+C IR*URR-Cl  I  *UR I 
GK  I  =  A 1 1  *UL R> A  1R*UL  I*BI*UCR*BR*UCUC1I*URR*C1R*URI 
RETURN 
ENC 


48 


II.  LISTING  OF  FORTRAN  IV  SUBROUTINE  FOR  ONE  STEP  FFT  ADVANCE 


SU6R0UT I NE  AOVANCI DRI*rv** i ah. tuLi/i 
DIMENSION  DR  II 2. 1 XM. JVM) *NN(2 1 »PV2( IXM) » WORK ( 130) 
DIMENSION  X I A( 64 ) f XRAI 64)  .  .  t 

REAL  PART  OF  D  IS  STORED  IN  DRH1.I.J) 

IMAGINARY  PART  OF  D  IS  STORED  IN  0RI(2.I.J) 

PV 2  CONTAINS  THE  SQUARE  OF  THE  FREQUENCIES  AS 
DEFINED  IN  EQ.I5-12).  THE  ASSUMPTION  IS  THAT 
J(  IN  REJ>.RT)«IXM  ANO  MIN  REPORT)*JVM.  IXM-JYM, 
DX*OY •  HENCE  P  DEFINED  BY  EQ.I5-12)  EQUALS 
Q  DEFINED  BY  EQ.I5-13). 

DZ02«.5*0Z 

ZOLD-VALUE  01  l  AT  BEGINNING  OF  STEP- THE  ONLY 

USE  OF  THIS  ARGUMENT  IS  TO  DETERMINE  FIRST  CALL 

FROM  SUBSEQUENT  CALLS 

FOURT-BRENNERS  FFT  SUBROUTINE-REF.  28 

NN ( 1  )*  IXM 

NN 1 2 1 *JYM 

CALL  F OURTI 0Rl*NN*2*-lf 1# WORK ) 

CHECK  F  >K  SKIPPING  OF  EVALUATION  OF  COMPLEX  ARRAY 

IF (ZOLO. EQ. .0 IGOTO  25 

IF ( ABSI  DZ02-DZ02‘J  ).GT..lE-7)G0Tu  25 

GOT  75 

00  c0  1=  . » IXM 

PH  I  =  -OZ  J2*PV2l  I) 

XRAI  I  |=C  )SI PHI ) 

X  I  A (  I )»S INC  PHI > 

CONTINUE 

00  100  1  =  1.  IXM 

CO  100  J=1.JYM 

OEN=  FLOATING  l)*NN<2)) 

YR*OR  III. I. J I/OEN 
Y I *0R I (  2. I» J l/DEN 
XR*XRAf I )*XKA(J)-X(A(I)*XIA(J) 

X  I  *XRA(  I  )*XIA(  JUXRAf  JI*XIA(  I  ) 

CRIIlf It J )»XR*YR-XIAVI 
0RI(2,I,  J|»XI*YRFXK*YI 
CONTINUE 

CALL  F  >UKT1  DR  I  »NNf  1.  1«  WORK ) 

DZ020«DZ  2 

RETURN 

ENC 


rJVM.DZ02.Z0LD) 


ES  AS 
THAT 
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LIST  OF  SYMBOLS 

B  Transformed  complex  wave  function 

Cp  Specific  heat  of  propagation  medium  at  constant  pressure 
D  Transformed  complex  wave  function 
F(£,n)  Spatial  distribution  of  beam  amplitude  at  aperture 
k  Wave  number  (=  2tt/A) 

L  Characteristic  length  in  propagation  direction  (=  R,  focused 
beams;  =  kr,  collimated  beams) 

n  Index  of  refraction  after  laser  heating  (function  of  £,  n,  and  c) 

nQ  Index  of  refraction  under  ambient  conditions  (assumed  constant) 

P  Total  beam  power 

r  Characteristic  length  in  £-n  plane;  e"1  folding  length  for  a 
Gaussian  beam 

R  Focal  distance  for  focused  beam 

T  Atmospheric  temperature  after  laser  heating  (function  of  £,  n, 
and  c) 

T0  Ambient  atmospheric  temperature 

U  Dimensionless  complex  wave  function  (=  W/Wq) 

V  Wind  speed  (assumed  to  be  in  the  positive  £  direction  and  may  be 
^-dependent) 

W  Complex  wave  function 

WQ  Normalizing  factor  [■  (P/nr2)*5] 

x,y  Adaptive  coordinates  in  transverse  plane  [=  X/a(Z),  Y/a(Z)] 

X,Y  Dimensionless  spatial  coordinates  in  transverse  plane  (=  5/r, 
n/r) 

Z  Dimensionless  coordinate  in  propagation  direction  (=  c/L) 
z  Transformed  coordinate  in  propagation  direction  [see  Eq.  (3-3)] 
aj  Absorption  coefficient 

a  Dimensionless  absorption  coefficient  (=ajL) 

6  Dimensionless  parameter  (=  kr2/L),  equal  to  the  Fresnel  number 
for  focused  beams  or  unity  for  collimated  beams 

Y  Dimensionless  thermal  blooming  distortion  parameter 

[=  (n02-l)P/(T0P0cpVAr)];  in  the  atmosphere  y  may  be  c-dependent. 

C  Spatial  coordinate  in  propagation  direction 

£,n  Spatial  coordinates  in  transverse  plane 

A  Laser  wavelength 
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p 

p, 
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Dimensionless  quantity  [*  j  kL(n2  -  nQ2)] 

Atmospheric  density  after  laser  heating  (function  of  £,  n,  s) 

Ambient  atmospheric  density 

Spatial  distribution  of  beam  phase  at  aperture 

Angular  slewing  velocity 

Dimensionless  slewing  velocity  (=  fiL/V) 
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